% clears unneceesary vriables after
% invoking readrts reads the header
% compute cal. table according to gain settings etc;
% 11.7.01 computes cross and auto spectra & average it 
% MT app rest & phse are now computed

% latest date 05.07.2002

clear;
readrts;

Notchi = NotchF(2:5);		% NOTE ! there may be problems 
IsNotch = Filt(2:5);		%        if notch do not have right dimension..
hpb1 = Filt(1);			% NOTE ! high pass filter is always disabled
prob = [Exd Eyd Eyd Eyd Eyd];	% 	 in calib05.m
gain = [Exa Eya Hxa Hya Hza];
clear Bden Record3 Bdend Bddm Count Record4 hpf Bdst Record5;
clear Bdstart MaxRecord ans nt Filt NBlocks array nlines Events ;
clear NotchF dummy Record1 factor Channel;

block = 1024;

notch = Notchi(find(IsNotch>0)); 
PAI = 3.1415926535897;

for i = 1:nchannel,
[cplx,frq1] = calib05(i,nband,block,prob(i),notch,hpb1,gain(i));
CalData(i,:) = cplx;
end;

%------------native declarations---------------

WI = hanning(block);
len = length(matrix)/block;
start = 100;
end1 = 160; 
st = start*block + 1;
en = (start+1)*block;
nchannel = 5;

[TLine,TLRad,TLFreq,TBand] = mtfreq(block); % target freq, line & radius for all bands

TLine1  = TLine(find(TBand==nband));
TLRad1  = TLRad(find(TBand==nband));
TLFreq1 = TLFreq(find(TBand==nband));

HxHx = zeros([length(TLine1),1]);
HyHy = HxHx;
HxHy = HxHx;
HyHx = HxHx;
ExHx = HxHx;
EyHx = HxHx;
ExHy = HxHx;
EyHy = HxHx;


for i = (start+1):end1,

%FFT, calibration etc.

for u = 1:5,
	K = CorTre(matrix(st:en,u),0);
	Channel1(:,u) = fft(K.*WI);
%Channel(:,u) = Channel1(1:block/2,u).*CalData(u,:)';
	Channel(:,u) = Channel1(TLine1,u).*conj(CalData(u,TLine1)');
end; 	



% auto and cross spectra---------------------------

HxHx = HxHx + Channel(:,3).*conj(Channel(:,3));
HyHy = HyHy + Channel(:,4).*conj(Channel(:,4));
HxHy = HxHy + Channel(:,3).*conj(Channel(:,4));
HyHx = HyHx + Channel(:,4).*conj(Channel(:,3));
ExHx = ExHx + Channel(:,1).*conj(Channel(:,3));
EyHx = EyHx + Channel(:,2).*conj(Channel(:,3));
ExHy = ExHy + Channel(:,1).*conj(Channel(:,4));
EyHy = EyHy + Channel(:,2).*conj(Channel(:,4));

%------------------------------------------------

clear Channel;

st = en + 1;
en = (i+1)*block;
end;




Zxy = HxHx.*ExHy - HxHy.*ExHx;
Zyx = HyHy.*EyHx - HyHx.*EyHy;
Den = HxHx.*HyHy - HxHy.*HyHx;


nfrq = length(TLine1);
frq = frq1;


for i = 1:nfrq,
   if abs(Den(i)>0)&abs(Zxy(i))>0&abs(Zyx(i))>0,
      ZXY(i)=Zxy(i)/Den(i);
      ZYX(i)=Zyx(i)/Den(i);
    end;
end;

for i = 1:nfrq,
   Phxy(i) = abs(atan(imag(ZXY(i))/real(ZXY(i))) * 180 / 3.1415926535897);
   Rhxy(i) = 0.2 * (1/TLFreq1(i)) * (ZXY(i) * conj(ZXY(i)));
   Phyx(i) = abs(atan(imag(ZYX(i))/real(ZYX(i))) * 180 / 3.1415926535897);
   Rhyx(i) = 0.2 * (1/TLFreq1(i)) * (ZYX(i) * conj(ZYX(i)));
end;


%--------------------------------------------------------------
T = 1./TLFreq1;
subplot(2,1,1);
hold on
loglog(T,Rhxy,'r*',T,Rhyx,'b*');
ylabel('ohm-m');
xlabel('period (s)');
axis([1e-2 1e2 1e-2 1e3]);
grid;
subplot(2,1,2);
hold on;
semilogx(T,Phxy,'r*',T,Phyx,'b*');
xlabel('period (s)');
ylabel('degrees');
title(file);
axis([1e-2 1e2 0 90]);

grid;

